Multi-parental fungal mapping population study to detect genomic regions associated with Pyrenophora teres f. teres virulence

In recent years multi-parental mapping populations (MPPs) have been widely adopted in many crops to detect quantitative trait loci (QTLs) as this method can compensate for the limitations of QTL analyses using bi-parental mapping populations. Here we report the first multi-parental nested association mapping (MP-NAM) population study used to detect genomic regions associated with host-pathogenic interactions. MP-NAM QTL analyses were conducted on 399 Pyrenophora teres f. teres individuals using biallelic, cross-specific and parental QTL effect models. A bi-parental QTL mapping study was also conducted to compare the power of QTL detection between bi-parental and MP-NAM populations. Using MP-NAM with 399 individuals detected a maximum of eight QTLs with a single QTL effect model whilst only a maximum of five QTLs were detected with an individual bi-parental mapping population of 100 individuals. When reducing the number of isolates in the MP-NAM to 200 individuals the number of QTLs detected remained the same for the MP-NAM population. This study confirms that MPPs such as MP-NAM populations can be successfully used in detecting QTLs in haploid fungal pathogens and that the power of QTL detection with MPPs is greater than with bi-parental mapping populations.

To understand the evolution of plant fungal pathogenicity, it is crucial to identify the genomic regions underlying the pathogenicity 1 . Initially, it was suggested that the pathogenicity of plant fungal pathogens was due to a single or a few genomic regions/genes and that these genes were recognised by the host plant and hence, follow a gene-for-gene model 2 . With the identification of several other host-pathogen interaction model systems, the expression of virulence of plant fungal pathogens on host plants was proposed to be a cumulative effect of multiple interacting genomic loci or quantitative trait loci (QTLs) 3,4 .
The QTLs responsible for inducing a virulence response in the host plant can be dissected through QTL mapping to extract important information on the genetic architecture of these phenotypic expressions 5,6 . QTL mapping studies are performed by linkage analysis in segregating bi-parental mapping populations 5 . The number of recombination events that occur within the mapping population determines the density of the linkage map. Many bi-parental mapping population studies have been conducted in fungal plant pathogens to identify QTLs responsible for the virulence of the pathogen [7][8][9][10] . However, in bi-parental mapping populations, the number of recombination events is restricted due to inadequate genetic diversity present between the two parents resulting in detecting broad chromosomal regions responsible for the phenotypic expression [11][12][13] . Hence, QTL mapping in bi-parental mapping populations is limited by low genetic map resolutions, inadequate genetic diversity between parents 14 and inability to identify the same QTL in other genetic backgrounds 15 .
To overcome the constraints of QTL mapping in bi-parental populations, genome-wide association mapping studies (GWAS) have been used to identify QTLs responsible for the virulence of fungal plant pathogens using genetically diverse populations from unknown kinship backgrounds 16 . The high mapping resolution in GWAS is achieved by its capacity to explore a broad genetic background coupled with capturing more allelic diversity occurring through historical recombination events and employing diversity panels which may comprise of fungal isolates collected from different geographical areas, years and hosts 11,16,17 . GWAS does, however, have its own Genotyping by DArTseq. Progeny  Phenotyping of the progeny isolates. Phenotyping of the progeny isolates was performed following a completely randomised design in a controlled environment room at the University of Southern Queensland, Australia, with three replicates, as described in Dahanayaka et al. 31 . The four barley cultivars (Beecher, Commander, Prior, and Skiff) were grown in pots with 5 cm diameter and 14 cm height with each pot containing four plants each of the four barley cultivars. Plants were grown at day (12 h) and night (12 h) temperature of 23 ± 1 °C and 17 ± 1 °C respectively, at 75% humidity for 14 days. Barley cultivars Beecher and Commander were used as the resistant and susceptible control, respectively.
The conidial suspensions for inoculations were prepared as described in Dahanayaka et al. 31 . Three millilitres of the suspension diluted to 10,000 conidia/mL were used per pot for inoculations 14 days after sowing. The parental isolate HRS10136 was used as the control isolate for each inoculation experiment to monitor disease reaction score differences across different experiments. After inoculation pots were incubated in the dark for 24 h at 95% humidity with a temperature of 20 ± 1 °C. Plants were then transferred to the controlled environment room with the same environmental conditions as mentioned above. Nine days after inoculation, disease reaction Genetic map construction. Individual genetic maps for four populations were constructed using Sili-coDArT and SNP marker data obtained from DArTseq™. Both markers were filtered using 10% as the cut-off value for the minimum amount of missing data per isolates. Non-polymorphic markers were removed along with markers deviating from the 1:1 segregation ratio. Clonal isolate pairs of each progeny were identified using the clonecorrect function in poppr package version 2.8.3 34 in R version 3.0.2 35 and one clonal isolate from each pair was removed. SilicoDArT and SNP markers were grouped into linkage groups using the make linkage groups function in MapManager QTXb20 version 2.0 36 with a p = 0.05 search linkage criterion. Markers were ordered using RECORD 37 and the final genetic map of each population was obtained after manual map curation 38 . The marker order of the linkage groups was confirmed by aligning marker positions to the Ptt reference genome, W1-1 (BioSample SAMEA4560035, BioProject PRJEB18107) using the bowtie2 39  QTLs associated with virulence in Ptt using the MP-NAM population. The genotypic and phenotypic data for the four bi-parental populations were combined to represent a single MP-NAM population. The final data set was quality filtered with 5% as the threshold minimum allele frequency for markers and 25% as the maximum missing data per isolate. This resulted in 1135 SilicoDArT and SNP markers from DArTseq™ and 399 progeny isolates being retained for the MP-NAM QTL analysis. QTL detection was carried out using cross-specific, bi-allelic and parental QTL effect models available in the mppR package. Significant thresholds for each trait were calculated using 1000 permutations. Cofactors were selected using the simple interval mapping method and multi-QTL model searches were implemented by composite interval mapping to detect QTLs in different types of effect models; bi-allelic (two alleles at the QTL position with consistent effect through the QTL position), cross-specific (at the QTL the allelic effects can be different in every cross), and parental (every parent carries a different allele with a consistent effect in every cross). The regression coefficients (R 2 ) for each QTL were detected. Cross validations for the detected QTLs were carried out to assess the QTL effect in a pseudoindependent population to confirm the putative QTLs.

Comparison between bi-parental and MP-NAM populations for QTL detection. The power of
QTL detection of bi-parental versus MP-NAM populations was compared by determining the number of QTLs identified with each set. Thus, results produced by the three MP-NAM QTL models were compared with results obtained from the four individual bi-parental mapping populations. Fifty isolates from each bi-parental population were randomly selected and pooled to create a smaller bi-parental and MP-NAM populations and the QTL analysis was conducted with Windows QTL Cartographer and the mppR package, respectively, following the same criteria as above. The analyses were repeated three times using three different sub-datasets (Run_1, Run_2 and Run_3). Each sub-dataset contained 50 randomly selected isolates from each bi-parental population.
Genetic diversity of the population. Principal analysis of coordinates (PCoA) and a neighbor-net network were used to detect the overall diversity structure and the genetic distribution of the isolates of each population. PCoA for the progeny isolates used in the MP-NAM population was conducted in DARwin v.6.0 45 using Euclidean distances with five principal coordinates. Neighbor-net network was built for the MP-NAM population in SplitsTree version 4.13 46 with 1000 bootstrap. The Neighbor-net network was built based on the method described by Bryant and Moulton 47 using neighbor joining algorithm 48 .

Plant ethical approval. Experimental research and field studies on plants (either cultivated or wild),
including the collection of plant material, must comply with relevant institutional, national, and international guidelines and legislation.

Results
DArTseq analysis. A total of 4809 SNPs and 9810 SilicoDArT were obtained from DArTseq™. After quality filtering of markers for 10% missing values, non-polymorphism and segregation distortion (1:1), 1241, 1102, 568 and 1168 SNPs and SlicoDArT markers were retained for the HRS11093xHRS09127, HRS11093xHRS10136, NB81xHRS09127 and NB81xHRS10136 populations, respectively and used for the phenotypic evaluation and genetic map construction ( Table 1).
Phenotyping of the progeny isolates. Disease reaction scores of the progeny isolates of the populations on barley cultivars Prior and Skiff ranged from avirulent to virulent ( Fig. 1 and Table 2). Some of the isolates used in the genetic map construction did not produce conidia for phenotyping, hence the difference in numbers. Disease reaction scores of these progeny isolates were skewed on the resistant control Beecher and the suscep- Genetic map and bi-parental QTL analyses. The genetic maps of the four bi-parental populations consisted of 12-16 linkage groups spanning from 2029 to 2683 cM ( Table 3). The average distance between flanking markers for the four populations ranged from 4.56 to 6.29 cM. The physical distance to genetic map distance ratio for the four populations with respect to the Ptt reference genome W1-1 ranged from 19.29 to 25.51 kb/cM. Results of the bi-parental mapping population QTL analysis are presented in Table 4 and Fig. 2. Between one to three QTLs associated with Prior or Skiff virulence were detected with individual bi-parental populations. A highly significant QTL associated with the Prior virulence was identified on chromosome 5 in all four populations. Significant LOD values for the QTLs ranged from 6.9 to 20.0 across the different populations. The phenotypic variance explained by the QTLs ranged from 33 to 63%. The most significant QTL identified for the Skiff virulence was located on chromosome 3 with a LOD score of 11.0 and 27% of the phenotypic variance explained. This QTL was only identified in the NB81xHRS09127 population. A QTL associated with Skiff virulence was detected on chromosome 5 in the NB81xHRS10136 population with a LOD score of 6.4 and 33% of the phenotypic variation explained. This QTL was located in the same region as the Prior virulence QTL.  cross-specific and parental QTL effect models using 399 progeny isolates. Maximum of five and three QTLs were identified for the Prior and Skiff virulence, respectively. All models identified a QTL on chromosome 5 for Prior virulence with LOD scores ranging from 30.9 to 36.7 ( Table 5). The phenotypic variance explained by the QTLs ranged from 14 to 39%. With both the cross-specific and parental models the QTL on chromosome 5 was split into two different QTLs. For Skiff virulence, three QTLs were detected across three different chromosomes with LOD scores ranging from 3.8 to 9.9 and phenotypic variance explained by these QTLs ranged from 4 to 13% across the different models. The same QTL was detected on chromosome 3 across all three models.  (Table 5). The QTLs detected with all three models for Prior virulence on chromosomes 4 and 5 in the original MP-NAM mapping population were also reported for all runs and models in the smaller MP-NAM population. Similarly, the chromosome 3 QTL associated with Skiff virulence was detected with all models and runs. Minor QTLs were only detected in some runs or not at all with some of the models, e.g. QTL on chromosome 8 associated with Prior virulence (Table 5). Using the bi-parental mapping populations only the QTL on chromosome 5 of the   www.nature.com/scientificreports/ HRS11093xHRS09127 population was detected in all runs ( Table 5). None of the QTLs associated with Skiff virulence were detected in all runs.
Genetic diversity of the population. The PCoA of 399 progeny isolates showed clear clustering of the four bi-parental populations as expected. The first and second principal coordinate axes (PCoA1 and PCoA2) explained 12 and 6% of the variance, respectively, and separated clusters by population. Out of the four populations, NB81xHRS09127 showed the highest genetic diversity within the population while HRS11093xHRS10136 showed the lowest (Fig. 3). The neighbor-net network consisted of two main groups. One group contained progeny isolates from HRS11093xHRS10136 and NB81xHRS10136 while the other group contained progeny isolates from HRS11093xHRS09127 and NB81xHRS09127. The neighbor-net network constructed for the MP-NAM population was highly reticulated as expected (Fig. 4).

Discussion
Multi-parental mapping populations have been used for many crops to detect QTLs as they have been useful in providing insight into the genetic architecture of complex phenotypes like grain yield, plant height and biotic and abiotic resistance amongst other phenotypes 44 . However, as per the authors' knowledge usage of MPPs to identify QTLs in fungi has not previously been reported. This the first study to use a MP-NAM population to identify QTLs in a fungus. Bi-parental mapping populations have been extensively used in the detection of QTLs in crops 14,19 . However, there are limitations with bi-parental mapping populations compared to MPPs in the successful detection of significant QTLs 11,16 . The principal limitations of bi-parental mapping populations are low genetic diversity resulting from genetic bottlenecking due to the choice of two parents and limited events of effective recombination to develop a precise genetic map. These limitations reduce the number of QTLs captured by bi-parental mapping populations.
The majority of NAM populations have been developed by crossing one reference parent with several other donor parents to obtain a diverse genetic background 49 . The NAM population used in this study was developed using more than one reference parent. The genetic principle behind the MP-NAM analysis used in this study is based on bi-parental populations and diversity panels. Bi-parental population analysis has the confounding effect of population structure, which is avoided when using MP-NAM analysis. Another advantage of MP-NAM populations is the detection of multiple alleles in the same QTL along with rarer and more diverse alleles.
Two out of the four parental isolates used in our mapping populations demonstrated high and low virulence on barley cultivars Prior and Skiff, respectively, while the other two demonstrated low and high virulence on Prior and Skiff, respectively. An increased number of recombination events occurring in a fungal population increases the genetic diversity of the population 50 . Using more than one parental isolate for the desired phenotype increases genetic diversity of the population as shown in the neighbour-net network and PCoA analysis in the current study and captures QTLs associated with the desired phenotype from more than one parent. As shown in this study, one to three QTLs were detected that were associated with Prior virulence using bi-parental mapping populations while five QTLs were detected using the MP-NAM population. For Skiff virulence, only two QTLs were reported for each bi-parental mapping population while three QTLs were detected for the MP-NAM population, thus confirming that the MP-NAM QTL detection method can detect a greater number of QTLs than the bi-parental method.
Most QTLs identified in the four bi-parental mapping populations of this study were co-localised, e.g. QTLs on chromosomes 4 and 5 identified for the Prior virulence. In populations HRS11093xHRS09127 and HRS11093xHRS10136 the contributing parent for the QTL on chromosome 5 was HSR11093 and in NB81xHRS09127 and NB81xHRS10136, the contributing parent for the QTL was NB81 suggesting that these isolates possess the same genomic regions on chromosome 5 associated with Prior virulence. QTLs associated with Prior virulence identified on chromosome 1, 3, and 8 of population NB81xHRS09127, NB81xHRS10136 and HRS11093xHRS10136, respectively, were unique QTL, not detected in the other populations. The QTL on chromosome 8 associated with Skiff virulence of bi-parental mapping population was also co-located with the QTL on chromosome 8 for Prior virulence. The QTL associated with Skiff virulence on chromosome 5 also colocated with the QTL associated with Prior virulence on chromosome 5. Co-localization of QTL associated with Prior and Skiff virulence could be due to either the presence of tightly linked genes in the same genomic region or lack of effective events of recombination between the two regions to identify the two genes as separate QTLs.
The QTLs identified on chromosomes 1, 4 and 5 for Prior virulence from the multi-parent population based on biallelic, cross-specific and parental models were located at the same genomic regions. The QTL detected on chromosome 5 for Prior virulence based on the biallelic model was a broad QTL. This QTL was detected as two QTLs in the cross-specific and parental models suggesting that there are two QTLs in this genomic region. This also suggests that the cross-specific and parental models may have better resolution power than the biallelic model. All QTLs associated with Prior virulence and QTL on chromosome 8 for Skiff virulence of the MP-NAM population were co-located with the QTLs identified by the four bi-parental mapping populations. Co-localisation of these QTLs validates the robustness of the QTLs identified in our study.
QTLs identified in our study were also co-localised with QTLs reported in previous bi-parental and GWAS studies. The location of the QTL identified on chromosome 5 for Prior virulence was similar to the QTLs AvrHar 24 , PttBee2 (Koladia et al. 28 ), PttBee_5, PttSki_5, QTL11, QTL12 (Martin et al. 4 ) and USQV5 (Dahanayaka et al. 31 ) identified previously through bi-parental and GWAS studies. The QTL, PttBee2 28 and PttBee_5 4,29 associated with Beecher virulence were co-localised with our QTL on chromosome 5 associated with Prior virulence. The four parental isolates used in the current study were avirulent on barley cultivar Beecher and hence, used as the resistant control. Co-localisation of the Prior virulence QTL and the Beecher virulence QTL suggests that this genomic region codes for a gene/genes which would induce multiple reactions in different Ptt isolates. This may be due to the presence of different alleles of the same gene in different Ptt isolates or to possessing more than one virulence gene in the same genomic region. The QTL associated with Prior virulence on chromosomes 8 in the current study shared the same genomic region as PttTif2 28 . Three QTLs responsible for Skiff virulence detected in the current study on chromosomes 3, 7 and 10 were co-located with QTL QTL7, PttPri_7 4,29 and VR2 27 , respectively.
Comparison of the bi-parental and MP-NAM analyses revealed that MP-NAM QTL analysis had greater power to detect QTLs than bi-parental QTL mapping, even with a low number of progeny isolates. The number of progeny obtained from crossing fungal isolates can vary 51 . In some instances, crosses may produce insufficient numbers of ascospores for bi-parental mapping analyses. On such occasions, a MPP like MP-NAM is a useful alternative method for identifying QTLs in populations with low progeny numbers, where inter-connected multiple populations are available. www.nature.com/scientificreports/ In conclusion, this study reported the first MPP QTL analysis conducted with a fungal population. The biparental mapping QTL analyses detected one to three QTLs associated with Prior virulence and two with Skiff virulence for the four populations developed from crossing P. teres f. teres isolates. The MP-NAM population, developed by combining the bi-parental mapping populations, detected five and three QTLs for Prior and Skiff virulence, respectively. The comparison of bi-parental and MP-NAM analyses revealed that QTL detection power of a MP-NAM population was greater than that of a bi-parental fungal mapping population even with a low number of progeny.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on request.